10000 REM Calculation of vapor-liquid equilibrium for aqueous KCl solution
10050 REM Hovey et al. (1990)
10100 DEFDBL A-H, M-Z
10150 DIM HGKG(40),II(40),JJ(40),BP(10),BQ(10)
10200 DIM ATZ(4),ADZ(4),AAT(4),AAD(4)
10250 DIM QR(11),QT(10),QZR(9),QZT(9),V(10),A(8),C(18)
10300 DIM COEF(10,10),BBC(10),BC(10)
10350 GOSUB *BLOCKDATA
10400 H2O=18.01534:KCL=74.551#:R=8.3144#
10450 RHOC=.322#:VC=H2O/(RHOC*10)
10500 INPUT"Input temperature (deg C)";TT
10550 IF TT<300 THEN GOTO 10500
10600 T=TT+273.15#
10650 TTT=T
10700 RT=R*T
10750 INPUT "Input molality of KCl in liquid phase";M
10760 IF M=<0 THEN GOTO 10750
10800 LPRINT USING"Temperature(deg C)=+###.###     molality of KCl=+##.######";TT,M
10850 LPRINT
10900 WLKCL=100*KCL*M/(1000+KCL*M)
10950 S=6509050/T/T-29128.5#/T+41.2623#+.0759889*T
11000 PVLS=-308.19+5.01055*TT-.0313178*TT*TT+9.29982D-005*TT*TT*TT
11050 PVLS=PVLS-1.19918D-007*TT*TT*TT*TT+5.46864D-011*TT*TT*TT*TT*TT
11100 IF WLKCL<S THEN GOTO 11350
11150 LPRINT"Supersaturated (or saturated) with KCl"
11200 LPRINT USING"Pressure = +#.###^^^^^  Solubility=+##.##";PVLS,1000*S/((100-S)*KCL)
11250 LPRINT : LPRINT : LPRINT : LPRINT
11300 GOTO 10500
11350 XL=M/(M+1000/H2O)
11400 YL=M*H2O/1000
11450 INPUT"Do you input intial estimates for DRL, DRV, and YV?  Input Y (or y) or N (n)";A$
11500 IF A$="N" OR A$="n" THEN GOTO 11950
11550 INPUT"DRL";DRL:INPUT"DRV";DRV:INPUT"YV";YV
11600 DELDL=1D-006:DELDV=1D-006
11650 MYULS=0:MYUVS=0:MYULWATER=0:MYUVWATER=0
11700 REM Parameters of Hovey et al.
11750 B10=116989.96#-157.8381#*T+.06641785*T*T-3.039773*10000000/T
11800 B11=-4452.32#+2.11429#*T-1.958284#*100000000000000/(T*T*T*T)
11850 B20=-37956.21#+45.32167#*T-1.831384#*1D+020/(T*T*T*T*T*T)
11900 GOTO 14450
11950 TC=647.096
12000 IF TT=<350 THEN GOTO 12100
12050 TAU=1#-623.15#/TC:GOTO 12150
12100 TAU=1#-T/TC
12150 GOSUB *WPKCL
12200 DRL0=DAQ*(1#-WLKCL/100)/RHOC:DRV0=DWATERV/RHOC
12250 XV=1D-012:YV=XV/(1#+XV)
12300 DRL=0:DRV=0
12350 DELDL=1D-006:DELDV=1D-006
12400 MYULS=0:MYUVS=0:MYULWATER=0:MYUVWATER=0
12450 REM Parameters of Hovey et al.
12500 B10=116989.96#-157.8381#*T+.06641785*T*T-3.039773*10000000/T
12550 B11=-4452.32#+2.11429#*T-1.958284#*100000000000000/(T*T*T*T)
12600 B20=-37956.21#+45.32167#*T-1.831384#*1D+020/(T*T*T*T*T*T)
12650 DRL=DRL0
12700 DRV=DRV0
12750 MYULS=RT*LOG(XL)+VC*((-1#)*B10/DRL+B11*(LOG(DRL)+1#/DRL))-2#*VC*YL*B20/DRL
12800 MYUVS=RT*LOG(YV)-RT*LOG(1#+YV)+VC*((-1#)*B10/DRV+B11*(LOG(DRV)+1#/DRV))-2#*VC*YV*B20/DRV
12850 DIFMYUS=MYUVS-MYULS
12900 DMYUSDY=RT/YV/(1#+YV)-2#*VC*B20/DRV
12950 IF ABS(DIFMYUS)<1D-010 THEN GOTO 13250
13000 IF ABS(DIFMYUS/DMYUSDY)<1D-021 THEN GOTO 13250
13050 YVDEC=DIFMYUS/DMYUSDY
13100 IF YVDEC<YV THEN GOTO 13200
13150 YVDEC=YVDEC*.5#:GOTO 13100
13200 YV=YV-YVDEC:GOTO 12800
13250 XV=YV/(1#+YV)
13300 X=0:PRES=0:PL=0:PV=0:PLD=0
13350 RHOV=DRV*RHOC
13400 X=RHOV
13450 GOSUB *WATER
13500 PV=PRES
13550 PV=PV+YV*(B10+B11*(DRV-1#))+YV*YV*B20
13600 RHOL=DRL*RHOC
13650 X=RHOL
13700 GOSUB *WATER
13750 PL=PRES
13800 PL=PL+YL*(B10+B11*(DRL-1#))+YL*YL*B20
13850 DIFPLV=PL-PV
13900 RHOLD=(DRL+DELDL)*RHOC
13950 X=RHOLD
14000 GOSUB *WATER
14050 PLD=PRES
14100 PLD=PLD+YL*(B10+B11*((DRL+DELDL)-1#))+YL*YL*B20
14150 DPDDL=(PLD-PL)/DELDL
14200 IF ABS(DIFPLV/DPDDL)<1D-010 THEN GOTO 14450
14250 DRLDEC=DIFPLV/DPDDL
14300 IF DRL>DRLDEC THEN GOTO 14400
14350 DRLDEC=.5#*DRLDEC:GOTO 14300
14400 DRL=DRL-DRLDEC:GOTO 13600
14450 GOSUB *SOLUTION
14500 LPRINT
14550 DRL=BC(1):DRV=BC(2):YV=BC(3)
14600 GOSUB *FUNC
14650 RHOLIQUID=(H2O+YL*KCL)*RHOC*DRL/H2O
14700 RHOVAPOR=(H2O+YV*KCL)*RHOC*DRV/H2O
14750 LPRINT USING"Pressure  liquid=+#.###^^^^^  vapor=+#.###^^^^^";PL,PV
14800 LPRINT USING"Density   liquid=+#.###^^^^^  vapor=+#.###^^^^^";RHOLIQUID,RHOVAPOR
14850 LPRINT:LPRINT:LPRINT:LPRINT
14900 GOTO 10500
14950 *WPKCL
15000 REM Densities of liquid and vapor under vapor-saturation
15050 WL1=1.99274064:WL2=1.09965342:WL3=-.510839303:WL4=-1.75493479#:WL5=-45.5170352:WL6=-674694.45#
15100 WV1=-2.0315024:WV2=-2.6830294:WV3=-5.38626492#:WV4=-17.2991605:WV5=-44.7586581#:WV6=-63.9201063
15150 DL=1#+WL1*EXP((1#/3#)*LOG(TAU))+WL2*EXP((2#/3#)*LOG(TAU))+WL3*EXP((5#/3#)*LOG(TAU))+WL4*EXP((16#/3#)*LOG(TAU))+WL5*EXP((43#/3#)*LOG(TAU))+WL6*EXP((110/3#)*LOG(TAU))
15200 DL=RHOC*DL
15250 DV=WV1*EXP((2#/6#)*LOG(TAU))+WV2*EXP((4#/6#)*LOG(TAU))+WV3*EXP((8#/6#)*LOG(TAU))+WV4*EXP((18#/6#)*LOG(TAU))+WV5*EXP((37#/6#)*LOG(TAU))+WV6*EXP((71#/6#)*LOG(TAU))
15300 DV=RHOC*EXP(DV):DWATERV=DV
15350 REM Haas equation
15400 BA1=-940.409:BA2=1716.28#:BA3=-774.341#:BA4=-126.905:BA5=116.456#
15450 VCON=3.106
15500 V0=1#/DL
15550 D1=1000*V0+M*BA1+M*BA2*V0+M*BA3*V0*V0+(BA4+BA5*V0)*V0*V0*M*SQR(M)/((VCON-V0)*(VCON-V0))
15600 DAQ=(1000+KCL*M)/D1
15650 RETURN
15700 *WATER
15750 GIB=0
15800 RTG=GASCON*T
15850 REM BBT
15900 V(1)=1#
15950 FOR I=2 TO 10
16000 V(I)=V(I-1)*TZ/T
16050 NEXT I
16100 B1=BP(1)+BP(2)*LOG(1#/V(2))
16150 B2=BQ(1)
16200 FOR I=3 TO 10
16250 B1=B1+BP(I)*V(I-1)
16300 B2=B2+BQ(I)*V(I-1)
16350 NEXT I
16400 DD=X
16450 D=DD*FD
16500 REM QQTD
16550 QR(1)=0
16600 Q5=0
16650 Q=0
16700 AR=0
16750 E=EXP(-AA*D)
16800 Q10=D*D*E
16850 Q20=1#-E
16900 QR(2)=Q10
16950 V=TZ/T
17000 QT(1)=T/TZ
17050 FOR I=2 TO 10
17100 QR(I+1)=QR(I)*Q20
17150 QT(I)=QT(I-1)*V
17200 NEXT I
17250 FOR I=1 TO INC
17300 K=II(I)+1
17350 L=JJ(I)
17400 QK=CDBL(K)
17450 QZR(K-1)=QR(K+1):QZT(L)=QT(L+1):QZR(K)=QR(K+2):QZT(L+1)=QT(L+2)
17500 QP=HGKG(I)*AA*QZR(K-1)*QZT(L)
17550 Q=Q+QP
17600 AR=AR+HGKG(I)*QZR(K)*QZT(L)/(Q10*QK*RTG)
17650 NEXT I
17700 QP=0
17750 FOR J=37 TO 40
17800 IF HGKG(J)=0 THEN GOTO 18900
17850 K=II(J)
17900 KM=JJ(J)
17910 QK=CDBL(K) : QKM=CDBL(KM)
17950 DDZ=ADZ(J-36)
18000 DEL=D/DDZ-1#
18050 IF ABS(DEL)<1D-010 THEN DEL=1D-010
18150 EX1=-AAD(J-36)*(DEL^QK)
18200 DEX=EXP(EX1)*(DEL^QKM)
18250 ATT=AAT(J-36)
18300 TX=ATZ(J-36)
18350 TAU=T/TX-1#
18400 EX2=-ATT*TAU*TAU
18450 TEX=EXP(EX2)
18500 Q10=DEX*TEX
18550 QM=QKM/DEL-QK*AAD(J-36)*(DEL^(QK-1#))
18600 FCT=QM*D*D*Q10/DDZ
18650 QP=QP+HGKG(J)*FCT
18700 AR=AR+Q10*HGKG(J)/RTG
18750 NEXT J
18800 Q=Q+QP
18850 REM BASEDT
18900 Y=.25#*B1*D
18950 XX=1#-Y
19000 Z0=(1#+G1*Y+G2*Y*Y)/(XX*XX*XX)
19050 Z=Z0+4#*Y*(B2/B1-GF)
19100 AB=-LOG(XX)-(G2-1#)/XX+28.16666667#/(XX*XX)+4#*Y*(B2/B1-GF)+15.166666667#+LOG(D*RTG/.101325)
19200 BASEF=Z
19300 PRES=FP*(RTG*D*BASEF+Q)
19350 Z=BASEF+Q/(RTG*D)
19400 REM THERMDT
19450 TIDEAL=T/100
19500 TL=LOG(TIDEAL)
19550 GI=-(C(1)/TIDEAL+C(2))*TL
19600 FOR I=3 TO 18
19650 GI=GI-C(I)*(TIDEAL^CDBL((I-6)))
19700 NEXT I
19750 AI=GI-1#
19800 AD=AB+AR+AI-UREF/T+SREF
19850 GD=AD+Z
19900 GIB=GD*RTG*FH
19950 RETURN
20000 *FUNC
20050 FSUMI=0
20100 GIB=0:PRES=0:PV=0:PL=0:X=0:DIFP=0:DIFMYUS=0:DIFMYUWATER=0
20150 DRL=BBC(1):DRV=BBC(2):YV=BBC(3)
20200 RHOV=DRV*RHOC
20250 X=RHOV
20300 GOSUB *WATER
20350 GIBV=GIB
20400 PV=PRES
20450 PV=PV+YV*(B10+B11*(DRV-1#))+YV*YV*B20
20500 RHOL=DRL*RHOC
20550 X=RHOL
20600 GOSUB *WATER
20650 GIBL=GIB
20700 PL=PRES
20750 PL=PL+YL*(B10+B11*(DRL-1#))+YL*YL*B20
20800 DIFP=PV-PL
20850 RHOVD=(DRV+DELDV)*RHOC
20900 X=RHOVD
20950 GOSUB *WATER
21000 GIBVD=GIB
21050 PVD=PRES
21100 PVD=PVD+YV*(B10+B11*((DRV+DELDV)-1#))+YV*YV*B20
21150 DPDDV=(PVD-PV)/DELDV
21200 DPDY=B10+B11*(DRV-1#)+2#*YV*B20
21250 RHOLD=(DRL+DELDL)*RHOC
21300 X=RHOLD
21350 GOSUB *WATER
21400 GIBLD=GIB
21450 PLD=PRES
21500 PLD=PLD+YL*(B10+B11*((DRL+DELDL)-1#))+YL*YL*B20
21550 DPDDL=(PLD-PL)/DELDL
21600 MYULS=RT*LOG(XL)+VC*((-1#)*B10/DRL+B11*(LOG(DRL)+1#/DRL))-2#*VC*YL*B20/DRL
21650 MYUVS=RT*LOG(YV)-RT*LOG(1#+YV)+VC*((-1#)*B10/DRV+B11*(LOG(DRV)+1#/DRV))-2#*VC*YV*B20/DRV
21700 DIFMYUS=MYUVS-MYULS
21750 DMYUSDDL=VC*B11/DRL+VC*(B10-B11+2#*YL*B20)/(DRL*DRL)
21800 DMYUSDDV=VC*B11/DRV+VC*(B10-B11+2#*YV*B20)/(DRV*DRV)
21850 DMYUSDY=RT/YV/(1#+YV)-2#*VC*B20/DRV
21900 MYULWATER=GIBL+YL*VC*(B10/DRL+B11*(1#-1/DRL))+2#*VC*YL*YL*B20/DRL+RT*LOG(1#-XL)
21950 MYUVWATER=GIBV+YV*VC*(B10/DRV+B11*(1#-1/DRV))+2#*VC*YV*YV*B20/DRV-RT*LOG(1#+YV)
22000 DIFMYUWATER=MYUVWATER-MYULWATER
22050 DMYUWATERDDL=(GIBLD-GIBL)/DELDL-YL*VC*(B10-B11+2#*YL*B20)/(DRL*DRL)
22100 DMYUWATERDDV=(GIBVD-GIBV)/DELDV-YV*VC*(B10-B11+2#*YV*B20)/(DRV*DRV)
22150 DMYUWATERDY=VC*(B10/DRV+B11*(1#-1#/DRV))+4#*YV*VC*B20/DRV-RT/(1#+YV)
22200 FSUMI=(PV/PL-1#)*(PV/PL-1#)+(MYUVS/MYULS-1#)*(MYUVS/MYULS-1#)+(MYUVWATER/MYULWATER-1#)*(MYUVWATER/MYULWATER-1#)
22250 RETURN
22300 *SOLUTION
22350 REM Non-linear equation
22400 REM Main routine
22450 REM Unknowns=KK
22500 EPS=1D-005:IMAX=100:AL=1#:KK=3
22550 DRLMAX=3.5#:DRLMIN=1.5#:DRVMAX=1#:DRVMIN=.02:YVMIN=1D-012
22600 YVMAX=YL-.001
22650 IF TT=<350 THEN GOTO 22750
22700 DRLMIN=1#:DRVMAX=1.5#:GOTO 22750
22750 LPRINT USING"Eps=+#.##^^^^^";EPS
22800 LPRINT USING"DRLMAX=+#.##  DRLMIN=+#.##  DRVMAX=+#.##  DRVMIN=+#.##  YVMAX=+#.##^^^^^  YVMIN=+#.##^^^^^";DRLMAX,DRLMIN,DRVMAX,DRVMIN,YVMAX,YVMIN
22850 BC(1)=DRL:BC(2)=DRV:BC(3)=YV
22900 LPRINT"Initial estimate"
22950 FOR I=1 TO KK
23000 LPRINT USING"BC(#)=+#.####^^^^^  ";I,BC(I);
23050 NEXT I
23100 LPRINT
23150 FOR I=1 TO KK:BBC(I)=BC(I):NEXT I
23200 GOSUB *FUNC
23250 LPRINT USING"Initial FSUM=+#.######^^^^^";FSUMI
23300 LPRINT
23350 FSUM=FSUMI
23400 ICON=KK:ITER=1
23450 GOSUB *BCSOLV
23500 LPRINT USING"ICON=##  FSUM=+#.######^^^^^    ITERATION=###";ICON,FSUM,ITER
23550 ITER=ITER+1:AL=1#
23600 FOR J=1 TO KK:LPRINT USING"BC(#)= +#.####^^^^^   ";J,BC(J);:NEXT J:LPRINT
23650 IF ITER>IMAX THEN ICON=-1
23700 LPRINT
23750 IF ICON>0 THEN GOTO 23450
23800 IF ICON=0 THEN GOTO 23950
23850 IF ICON=-1 THEN GOTO 23900
23900 LPRINT "Function improvement possible"
23950 LPRINT "Solution"
24000 FOR J=1 TO KK
24050 LPRINT USING"BC(#)= +#.##########^^^^^";J,BC(J)
24100 NEXT J
24150 RETURN
24200 REM Subroutine BCSOLV
24250 *BCSOLV
24300 KP1=KK+1
24350 FOR I=1 TO KK:FOR J=1 TO KP1:COEF(I,J)=0:NEXT J:NEXT I
24400 FOR I=1 TO KK:BBC(I)=BC(I):NEXT I
24450 GOSUB *FUNC
24500 COEF(1,KP1)=(-1#)*DIFP:COEF(2,KP1)=(-1#)*DIFMYUS:COEF(3,KP1)=(-1#)*DIFMYUWATER
24550 COEF(1,1)=(-1#)*DPDDL:COEF(1,2)=DPDDV:COEF(1,3)=DPDY
24600 COEF(2,1)=(-1#)*DMYUSDDL:COEF(2,2)=DMYUSDDV:COEF(2,3)=DMYUSDY
24650 COEF(3,1)=(-1#)*DMYUWATERDDL:COEF(3,2)=DMYUWATERDDV:COEF(3,3)=DMYUWATERDY
24700 REM Solve correction equation
24750 FOR ISL1=1 TO KK
24800 ISL2=ISL1+1
24850 FOR ISL3=ISL2 TO KP1
24900 COEF(ISL1,ISL3)=COEF(ISL1,ISL3)/COEF(ISL1,ISL1)
24950 NEXT ISL3
25000 FOR ISL3=1 TO KK
25050 IF ISL1<>ISL3 THEN GOTO 25150
25100 GOTO 25300
25150 FOR ISL4=ISL2 TO KP1
25200 COEF(ISL3,ISL4)=COEF(ISL3,ISL4)-COEF(ISL1,ISL4)*COEF(ISL3,ISL1)
25250 NEXT ISL4
25300 NEXT ISL3
25350 NEXT ISL1
25400 FOR J1=1 TO KK
25450 BBC(J1)=BC(J1)+COEF(J1,KP1)
25500 NEXT J1
25550 ICON=0
25600 IF BBC(1)>DRLMIN THEN GOTO 25700
25650 BBC(1)=DRLMIN:GOTO 25800
25700 IF BBC(1)<DRLMAX THEN GOTO 25800
25750 BBC(1)=DRLMAX
25800 IF BBC(2)>DRVMIN THEN GOTO 25900
25850 BBC(2)=DRVMIN:GOTO 26000
25900 IF BBC(2)<DRVMAX THEN GOTO 26000
25950 BBC(2)=DRVMAX
26000 IF BBC(3)>YVMIN THEN GOTO 26100
26050 BBC(3)=YVMIN:GOTO 26200
26100 IF BBC(3)<YVMAX THEN GOTO 26200
26150 BBC(3)=YVMAX
26200 GOSUB *FUNC
26250 IF FSUMI<FSUM THEN FSUM=FSUMI:GOTO 27400
26300 AL=AL*1#/2#
26350 FOR J1=1 TO KK
26400 BBC(J1)=BC(J1)+AL*COEF(J1,KP1)
26450 NEXT J1
26500 IF BBC(1)>DRLMIN THEN GOTO 26600
26550 BBC(1)=DRLMIN:GOTO 26700
26600 IF BBC(1)<DRLMAX THEN GOTO 26700
26650 BBC(1)=DRLMAX
26700 IF BBC(2)>DRVMIN THEN GOTO 26800
26750 BBC(2)=DRVMIN:GOTO 26900
26800 IF BBC(2)<DRVMAX THEN GOTO 26900
26850 BBC(2)=DRVMAX
26900 IF BBC(3)>YVMIN THEN GOTO 27000
26950 BBC(3)=YVMIN:GOTO 27100
27000 IF BBC(3)<YVMAX THEN GOTO 27100
27050 BBC(3)=YVMAX
27100 GOSUB *FUNC
27150 IF FSUMI<FSUM THEN FSUM=FSUMI:GOTO 27400
27200 AL=AL*.5#
27250 IF AL>1D-008 OR AL<-1D-008 THEN GOTO 26350
27300 IF AL>=-1D-008 AND AL<0 THEN LPRINT "Failed. ":END
27350 AL=-1:GOTO 26350
27400 FOR I=1 TO KK
27450 IF ABS(COEF(I,KP1)/BBC(I))>=EPS THEN ICON=ICON+1
27500 NEXT I
27550 FOR I=1 TO KK:BC(I)=BBC(I):NEXT I
27600 RETURN
27650 *BLOCKDATA
27700 FOR I=1 TO 4:READ ATZ(I):NEXT I
27750 DATA 640#,640#,641.6#,270#
27800 FOR I=1 TO 4:READ ADZ(I):NEXT I
27850 DATA 0.319#,0.319#,0.319#,1.55#
27900 FOR I=1 TO 4:READ AAT(I):NEXT I
27950 DATA 2.0D+004,2.0D+004,4.0D+004,25.0#
28000 FOR I=1 TO 4:READ AAD(I):NEXT I
28050 DATA 34.0#,40.0#,30.0#,1.050#
28100 GASCON=.461522#:TZ=647.073:AA=1#:INC=36
28150 UREF=-4328.454977#:SREF=7.6180720#
28200 G1=11#:G2=44.333333333333#:GF=3.5#
28250 FOR I=1 TO 10:READ BP(I):NEXT I
28300 DATA 0.7478629#,-0.3540782#,0.0#,0.0#,0.007159876#,0.0#,-0.003528426#,0.0#,0.0#,0.0#
28350 FOR I=1 TO 10:READ BQ(I):NEXT I
28400 DATA 1.1278334#,0.0#,-0.5944001#,-5.010996#,0.0#,0.63684256#,0.0#,0.0#,0.0#,0.0#
28450 FOR I=1 TO 40:READ HGKG(I):NEXT I
28500 DATA -5.3062968529023D+002,2.2744901424408D+003,7.8779333020687D+002
28550 DATA -6.9830527374994D+001,1.7863832875422D+004,-3.9514731563338D+004
28600 DATA 3.3803884280753D+004,-1.3855050202703D+004,-2.5637436613260D+005
28650 DATA 4.8212575981415D+005,-3.4183016969660D+005,1.2223156417448D+005
28700 DATA 1.1797433655832D+006,-2.1734810110373D+006,1.0829952168620D+006
28750 DATA -2.5441998064049D+005,-3.1377774947767D+006,5.2911910757704D+006
28800 DATA -1.3802577177877D+006,-2.5109914369001D+005,4.6561826115608D+006
28850 DATA -7.2752773275387D+006,4.1774246148294D+005,1.4016358244614D+006
28900 DATA -3.1555231392127D+006,4.7929666384584D+006,4.0912664781209D+005
28950 DATA -1.3626369388386D+006,6.9625220862664D+005,-1.0834900096447D+006
29000 DATA -2.2722827401688D+005,3.8365486000660D+005,6.8833257944332D+003
29050 DATA 2.1757245522644D+004,-2.6627944829770D+003,-7.0730418082074D+004
29100 DATA -0.225#,-1.68#,0.055#,-93.0#
29150 FOR I=1 TO 40:READ II(I):NEXT I
29200 DATA 0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6,8,8,8,8,2,2,0,4,2,2,2,4
29250 FOR I=1 TO 40:READ JJ(I):NEXT I
29300 DATA 2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,1,4,4,4,0,2,0,0
29350 FOR I=1 TO 18:READ C(I):NEXT I
29400 DATA 1.9730271018D+001,2.09662681977D+001,-4.83429455355D-001,6.05743189245D+000,2.256023885D+001,-9.87532442D+000,-4.3135538513D+000,4.58155781D-001
29450 DATA -4.7754901883D-002,4.1238460633D-003,-2.7929052852D-004
29500 DATA 1.4481695261D-005,-5.6473658748D-007,1.6200446D-008,-3.303822796D-010,4.51916067368D-012,-3.70734122708D-014,1.37546068238D-016
29550 FD=1#:FP=10:FH=18.01534
29600 RETURN